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Nonlinear  Optical  Fiber  Simulation 

A  paper  documenting  our  developments  in  three-dimensional  optical  fiber  simulation  has 
appeared  in  the  literature.  A  graduate  student,  Jun  Liu,  used  these  developments  to 
simulate  three-dimensional  soliton  simulation  in  an  optical  fiber.  This  was  his  master’s 
thesis. 

Quantum  Semiconductor  Simulation 

The  quantum  semiconductor  simulation  project  has  been  pursued  with  the  intention  of 
developing  quantum  switching  mechanisms.  This  is  being  done  in  collaboration  with 
the  Quantum  Optics  Theory  Group  of  the  Physics  Department  of  Washington  State 
University  headed  by  Professor  David  Citrin.  A  toe  domain  formulation  of  the 
Schroedinger  equation  has  been  developed.  The  specific  configuration  being  simulated 
is  a  quantum  dot.  This  was  chosen  broause  both  theoretical  and  experimental  data  are 
available  for  verification  of  the  accuracy  of  the  simulation  ( R.  C.  Ashoori,  “Electrons  in 
Artificial  Atoms,”  Nature,  379,  413-419,  1996).  Two  major  milestones  have  been 
achieved  so  far: 

Determination  of  the  Eigenstates  of  Arbitrary  Structures. 

The  analytical  description  of  particles  in  a  quantum  structure  is  available  only  for  the 
simplest  canonical  configurations.  Using  a  finite-difference  formulation  of  the 
Shroedinger  equations,  a  method  has  been  developed  to  determine  the  eigenenergies  and 
eigenfunctions  of  any  arbitrary  quantum  structure.  (Please  see  enclosed  manuscript.) 

Multiple  particles  in  a  quantum  dot 

The  interaction  of  two  electrons  in  a  quantum  dot  under  the  influence  of  a  magnetic  field 
has  been  simulated  using  Schroedinger  equation  and  the  Hartree-Folk  approximation. 
The  Hartree-Folk  formulation  takes  into  account  the  Coulomb  interaction  of  the  two 
particles  as  well  as  the  exchange  term,  which  is  a  purely  quantum  mechanical  effect 
based  on  the  relative  spins  of  the  two  particles.  The  Hartree-Folk  formulation  is 
computationally  very  intensive,  and  can  be  analytically  calculated  for  only  the  simplest 
cases.  By  reformulating  the  Coulomb  and  exchange  terms  as  convolutions,  a  fast  two- 
dimensional  Fourier  transform  algorithm  available  on  the  Cray  T90  of  the  San  Diego 
Supercomputer  Center  has  made  these  calculations  tractable.  Thus  far,  the  interaction  of 
two  particles  in  a  quantum  dot  subject  to  variable  external  magnetic  field  has  been 
simulated.  The  resulting  energy  levels  are  found  to  be  in  agreement  with  those  available 
in  the  literature.  (Please  see  the  enclosed  manuscript.) 
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Simulation  of  Terahertz  Pulse  Generation 

The  generation  of  very  short,  high  frequency  pulses  has  become  of  great  interest  to  the 
physics  community  over  the  past  few  years.  The  physical  processes  that  produce 
terahertz  (THz)  pulses  do  not  produce  the  short,  well  defined  pulses  for  applications 
such  as  time-domain  spectroscopy.  They  require  filtering.  One  such  filtering  technique 
is  spatiotemporal  shaping  in  which  the  pulse  is  passed  through  a  slot  in  a  metal  screen. 
This  is  an  application  well  suited  to  simulation,  where  different  apetoe  configurations 
can  be  tested  without  the  construction  of  expensive  and  time  consuming  experiments.  A 
paper  appeared  using  the  finite-difference  time-domain  method  to  do  a  two-dimensional 
simulation.  [J.  Bromage,  et  al,  “Spatiotemporal  shaping  of  half-cycle  terahertz  pulse  by 
diffraction  through  conductive  apertures  of  finite  thickness,”  /.  Optical  Soc.  Amer.  B, 
vol.  15,  pp  1399-1405,  April,  1998.]  Clearly  a  three-dirnensional  sirnulation  is  needed 
to  account  for  all  parameters.  However,  the  three-dimensional  simulation  of  the  far  field 
of  the  pulse  from  an  aperture  is  a  computationally  prohibitively  large  problem.  A 
method  was  developed  to  calculate  the  far  field  from  an  aperture  while  simulating  only 
the  immediate  area  around  the  aperture.  (This  is  described  in  a  paper.)  Furthermore,  the 
symmetry  of  the  problem  is  exploited  to  cut  the  computation  by  one  fourth.  Using  these 
methods  and  state-of-the-art  computer  resources,  the  three-dimensiond  simulation  of 
THz  pulse  was  accomplished.  A  graduate  student,  Sunil  Nekkanti,  has  used  this 
method  for  terahertz  pulse  shaping.  (Please  see  enclosed  manuscript.)  This  was  the 
subject  of  his  master’s  thesis. 

(5)  Technology  Transfer 

N/A 
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Abstract 

A  time-domain  simulation  method  is  presented  that  utilizes  the  Hartree-Fock  formulation 
to  characterize  two  particles  in  a  quantum  dot.  The  basis  of  the  simulation  is  the  finite- 
difference  time-domain  (FDTD)  method.  The  computation  is  made  tractable  by  formulating  the 
Coulomb  and  exchange  terms  as  digital  filtering  problems,  and  utilizing  two-dimensional  fast 
Fourier  transforms.  Two-electron  wavepacket  dynamics  are  calculated. 


This  material  is  based  upon  work  supported  by  the  U.  S.  Army  Research  Office  under  grant 
number  DAAH04-96- 1-0406,  and  by  a  grant  for  supercomputer  time  from  the  San  Diego 
Supercomputer  Center.  D.  S.  Citrin  was  supported  by  the  Office  of  Naval  Research  and  by  the 
National  Science  Foundation  through  grant  no.  DMR-9704503. 
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I.  Introduction 

Semiconductor  quantum  dots  (QD)  have  attracted  great  attention  due  to  the  quantization 
of  carriers  in  three  dimensions,  leading  to  discrete  spectra  [1].  Among  other  things,  they 
present  the  possibility  of  studying  the  detailed  interaction  of  particles  in  a  controlled 
environment  [2].  The  advent  of  measurement  techniques,  such  as  single-electron  capacitance 
spectroscopy  (SECS),  has  made  possible  the  determination  of  the  energy  of  individual  particles 
as  they  are  added  to  a  QD  under  various  conditions  [3, 4]. 

Along  with  the  growing  body  of  experimental  work,  there  have  been  efforts  to 
characterize  these  interactions  through  approximation  techniques  [5,  6,  7].  The  Hartree-Fock 
approximation  is  a  particularly  convenient  approach  for  the  interaction  of  multiple  particles  [6]. 
This  paper  presents  a  formulation  of  the  Hartree-Fock  approximation  using  the  finite-difference 
time-domain  (FDTD)  method.  The  FDTD  method  is  one  of  the  most  widely  used  methods  in 
electromagnetic  simulation  [8,  9]  and  it  has  recently  been  applied  to  the  simulation  of  the 
Schroedinger  equation  [9, 10].  In  this  paper,  the  simulation  of  two  electrons  is  presented.  This 
technique  allows  for  the  simulation  of  two-electron  wavepacket  dynamics  as  well  as  for  the 
determination  of  energy  eigenstates.  The  computational  requirements  of  the  Coulomb  and 
exchange  terms  are  partly  circumvented  by  using  signal-processing  techniques  and  a  two- 
dimensional  fast  Fourier  Transform  (FFT).  While  the  resulting  simulation  is  computationally 
intense,  it  is  well  within  the  realm  of  state-of-the-art  computing  platforms. 

Section  n  describes  the  FDTD  formulation  of  a  particle  in  a  two-dimensional  harmonic 
oscillator  in  a  magnetic  field.  This  is  the  usual  characterization  of  a  quantum  dot  [2].  Section 
in  shows  results  of  the  simulation  of  the  first  few  eigenstates  under  the  influence  of  a  magnetic 
field  and  verifies  their  accuracy  by  comparison  with  analytic  results.  In  Section  IV,  the 
implementation  of  the  Hartree-Fock  approximation  for  the  simulation  of  two  particles  is 
described  [11].  In  Section  V,  the  simulation  of  two  electrons  in  a  quantum  dot  is  presented. 

The  chemical  potential  of  the  first  two  electrons  is  found  to  be  in  excellent  agreement  with 
results  available  in  the  literature  [1,2]. 

II.  FDTD  Formulation  of  the  Schroedinger  Equation 
Basic  Formulation 

We  treat  a  QD  in  which  the  confinement  in  the  z  direction  is  much  stronger  than  in  the  x  and 
y  directions.  Thus,  the  energy-level  separation  associated  with  the  z  direction  quantization  is  much 
larger  than  any  other  energy  scale  in  the  problem.  We  therefore  assume  a  single  z-quantized  state. 
In  addition,  in  the  present  study  we  neglect  the  spin-orbit  interaction  which  is  relatively  small. 
(This  will  be  topic  of  a  future  study.)  The  time-domain  formulation  of  the  Schroedinger  equation 
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for  a  particle  in  a  two-dimensional  harmonic  oscillator  subject  to  a  magnetic  field  in  the 
perpendicular  direction  is  then  [12] 

=  — f- V -q-A]  yr(x,y)  +  V„o(x,y)-\i/(x,y),  (1) 

dt  2m  Vi  J 

where  V„o{,x,y)  is  the  two-dimensional  harmonic  oscillator  potential.  Equation  (1)  can  be 


written 


dy/ix,y)  _  ^ 

dt  2mh  [  dx^  dx^ 


i\ir{rx,y) 


VgV±Zl  +  2^1-4  +  x-f\w(x,y) 


dx  dy 


-iL^(x^+y‘r.f(x.y). 

it 

where  is  the  ground  state  energy  of  the  harmonic  oscillator  and  B  is  the  strength  of  the 
magnetic  field,  which  is  assumed  to  be  uniform  and  in  the  z  direction.  Next,  \j/  is  separated 
into  real  and  imaginary  parts: 

\l/(r,t)  =  WrealM  + 

Then  the  finite-difference  approximations  to  the  spatial  and  temporal  derivatives  are  taken. 
Equation  (2)  becomes  two  coupled  equations,  the  real  part  given  by 


yf  reality  j)  ~  realijij') 

h  At  imagii,  j)  yt"  /nKig(j  +  l>y)  “I"  imag(i~l,j) 

2m  Ax^  +  + 


+  ^^At\Ax^(i -  icf  +  Ax^U  -  jcfh"  ^'^imagiij) 
2mti  ^ 


(yf''~\eal(,Uj  +  1)  -  V^"  'realiij  “  1)) 

Ax(i  -  icy - — - 

2Ax 

+  1,;)  -  yf^^Keadi  “  1,7)) 

-  Ax(7  -  jc) - — 


(3) 


H - Uq  •  [(/  —  ic)  •  Ax  +  (7  —  7c)  •  Ax]  y/"  ^'^imag(i,j)- 

2r 

There  is  a  corresponding  equation  to  calculate  the  imaginary  part.  Here  we  discretize 
both  time  and  space;  the  time  step  is  indexed  by  n,  the  space  grid  by  i  and  7'.  Details  of  the 
above  derivation  are  given  in  [10]. 
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III.  Simulation  of  One  Particle  in  a  Quantum  Dot 

In  this  section,  we  will  use  the  techniques  described  in  the  previous  section  to  simulate 
the  lowest  order  eigenstates  in  a  quantum  dot.  Following  the  work  of  Ashoori  et  al.  [4],  we 
will  assume  the  dot  is  adequately  described  by  a  two-dimensional  harmonic  oscillator  with  a 
fundamental  energy  splitting  between  successive  single-particle  levels  of  5.4  meV.  The 
eigenstates  are  well  known  [12].  A  few  of  the  lower  order  states  are  shown  in  Fig.  1.  For 
simplicity,  we  only  plot  the  real  part  of  the  wavefunction  \|/.  We  describe  the  states  by  the 
quantum  numbers  (n,  1),  where  n,  the  principal  quantum  number,  is  a  positive  integer 
corresponding  to  the  number  of  nodes  in  the  wavefunction  moving  outward  from  the  center, 
and  I  is  the  axial  quantum  number  such  that  2  I/I  is  the  number  of  nodes  moving  in  a  circle  on  a 
constant  radius  around  the  dot  [2].  The  integer  I  can  be  positive  or  negative  corresponding  to  a 
waveform  that  is  moving  counterclockwise  or  clockwise,  respectively,  around  the  center  of  the 
dot.  The  energy  of  the  function  is  give  by  [2] 

E„,  =  nc0o(2n+\l\+l).  (4) 

Figure  2  is  a  time-lapsed  simulation  of  the  real  part  of  the  (0,1)  state.  Its  initial  energy  is 
10.8  meV,  as  calculated  by  Eq.  (4).  After  0.128  ps,  it  has  moved  one-third  of  a  revolution 
counterclockwise.  After  0.255  ps,  it  has  gone  two-thirds  of  a  revolution,  and  after  0.383  ps,  it 
has  returned  to  its  original  position.  This  is  not  surprising,  since  this  is  the  revival  time 
corresponding  to  10.8  meV.  This  simulation  used  cell  sizes  of  2  nm  and  time  steps  of  0.05  fs. 
The  total  simulation  space  was  60  by  60  cells  to  simulate  an  area  of  120  by  120  nm. 

Figure  3  is  a  similar  simulation  of  the  (0, 1)  state,  but  with  an  increasing  magnetic  field. 
Now  the  total  energy  is  calculated  by  the  formula  [2] 

+  +0)o-(2n+\l\+l),  (5) 

where  0)^  =  eBlm  is  the  cyclotron  frequency.  Note  that  when  the  particle  is  initialized  without 
the  magnetic  field,  the  energy  of  10.8  meV  is  evenly  divided  between  kinetic  and  potential 
energy.  The  radius  is  18.9  nm.  (The  values  of  energy  and  radius  are  actually  expectation 
values  of  the  corresponding  observables.)  As  the  magnetic  field  is  applied,  a  greater  portion  of 
the  energy  is  potential  energy,  and  the  radius  becomes  smaller.  These  two  phenomena 
correspond  to  the  additional  confinement  caused  by  the  magnetic  field.  The  total  energies 
correspond  closely  with  the  analytic  energy  calculated  by  Eq.(5). 

Such  a  simulation  was  made  for  each  of  the  six  lowest  energy  levels  and  is  plotted  in 
Fig.  4.  These  are  the  Fock-Darwin  levels  [1].  Results  are  shown  for  the  FDTD  simulation  and 
the  analytic  values  calculated  from  Eq.  (5).  In  general,  agreement  is  excellent.  Some 
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discrepancy  occurs  for  the  (0, 2)  level  at  magnetic  field  levels  above  6  T.  The  (0,  2)  function  is 
a  relatively  complex  one  (Fig.  1),  and  as  it  is  compressed  into  tighter  radii,  the  spatial  resolution 
of  2  nm  is  no  longer  adequate  to  maintain  accuracy.  Obviously,  this  could  be  overcome  by  a 
higher  resolution  and  a  larger  simulation  space. 

IV.  The  Hartree-Fock  Approximation 

The  Hartree-Fock  formulation  for  two  particles  results  in  the  following  coupled 
equations  [12]; 

ih  V,  -  « ■  a] V,  (/; )  +  V„  „.(/i  )i(^,  ('i ) 

dt  lm\i  J  ^ 

2  (6  a) 

-  /- .  WMS.,.,. 

|Ai-r2|  4;reo-'  r^-r^ 


=  if- V,  -  « ■  aIViC'-j)  +  Vn  aMViir,) 

*  2™'^-  J  (6b) 

^  ^  J,,  ^ ^  jA, 

d/ceoJ  ki-'il  4;r£o'’  ^-''21 

In  each  equation,  the  third  term  on  the  right  is  the  Coulomb  potential  and  the  last  term  is  the 
exchange  term.  is  the  Kroneker  delta  relating  the  spins  Sj  and  ^2  of  the  particles.  Only 

when  the  spins  are  in  the  same  direction  does  the  exchange  term  enter  into  the  calculations. 
Both  and  iff 2  are  complex  variables,  and  each  equation  results  in  two  separate  equations, 
similar  to  those  described  in  the  previous  section.  However,  they  are  coupled  by  the  Coulomb 
and  exchange  terms.  Note  that  each  of  these  integrals  is  a  spatial  integral  that  must  be  calculated 
for  each  position  in  the  problem  space.  This  threatens  to  overwhelm  the  computation.  This  is 
largely  overcome  by  recasting  these  integrals  as  convolutions  and  using  a  fast  two-dimensional 
Fourier  transform  for  the  calculation. 

We  begin  by  taking  the  Coulomb  integral  of  the  first  equation,  Eq.  (6a),  which  we  will 
call  /ci(q), 


^cM)  =  J dfz  =  J dr2  h(r^  -  r^)f{r^) 


where 


Kn  -r^)  = 


n-h\ 
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/(^2)  =  lV'2('‘2)f- 

This  can  be  written  in  the  xy  co-ordinates  as 

^ci(^i’3'i)  =  I ^2/ ^>’2  h{x\  -  x^,y^  -  yT)f{x2,y^) . 

Convolution  in  the  spatial  domain  gives  multiplication  in  the  k  domain,  so  taking  the  two- 
dimensional  Fourier  transforms  of  Eq.  (7)  gives 

kiK.K)  =  HKk,)F[K^,K)- 

We  have  eliminated  the  integration.  Of  course,  a  two-dimensional  Fourier  transform 
must  be  carried  out  at  every  time  step  to  get  because  the  corresponding  value  of 

fix2,y2)  =  lV'2('^2»>'2)f  is  updated  at  every  time  step;  then  an  inverse  transform  is  performed  to 
get  /(xpy,).  (  Hik^^  )  need  only  be  calculated  once  at  the  beginning  of  the  program.)  Fast 

Fourier  transforms  are  usually  implemented  in  specialized  subroutines  that  use  complex  arrays. 
Since  the  calculation  of  the  Coulomb  potential  involves  only  real  numbers,  both  Coulomb 
integrals  can  be  calculated  at  once,  one  in  the  real  buffer  and  one  in  the  imaginary  buffer. 

Now  we  turn  our  attention  to  the  calculation  of  the  exchange  terms.  The  two  exchange 
terms  of  Eqs.  (6  a)  and  (6  b)  are 

^Exl  ~  [ 


‘&2 


V'2*(^2¥i(^2) 

ki-^2l 

(8  a) 

Wx(ry)¥2iri) 

^1-^2 

(8  b) 

Note  that  rj  and  r,  in  Eqs.  (8  a)  and  (8  b)  are  just  integration  parameters.  If  they  were  both  set 
to  r,  for  instance,  then  we  can  see  that  (8  a)  and  (8  b)  are  just  complex  conjugates  of  each  other. 
Therefore,  it  is  only  necessary  to  calculate  one  integral  and  take  the  complex  conjugate  of  it  to 
get  the  other.  Of  course,  we  will  use  the  same  basic  approach  as  used  in  the  calculation  of  the 

Coulomb  potentials:  take  the  FFT  of  Y2ih)¥iih)^  multiply  it  by  the  transformed  version  of  the 

function  1  /  |r,  -  r2| ,  and  take  the  inverse  FFT  of  the  result. 

The  following  is  the  procedure  that  must  take  place  at  every  time  step: 

1)  One  forward  FFT  to  get  the  Fourier  transform  of  the  magnitudes  of  |vr,  (a;  )|^  and  |v'^2('2)r  ’  ^ 
multiplication  by  H{k^^  )>  then  an  inverse  FFT  to  get  the  two  Coulomb  integrals. 


2)  A  forward  FFT  of  multiplication  by  and  an  inverse  FFT  to  get 

the  exchange  term  for  The  exchange  term  for  1/^2  ('i)  is  just  the  complex  conjugate. 
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V.  Simulation  of  Two  Particles  in  a  Quantum  Dot 

In  this  section,  we  use  the  FDTD  formulation  of  the  Hartree-Fock  approximation 
described  in  the  previous  section  to  simulate  two  electrons  in  a  quantum  dot  as  described  by  the 
two-dimensional  harmonic  oscillator  potential.  We  also  use  the  same  configuration  in  the  FDTD 
program,  i.  e.,  a  cell  size  of  2  nm,  a  time  step  of  0.05  fs,  and  a  total  problem  space  that  is  60  by 
60  cells.  The  difference  is  that  there  are  two  simultaneous  simulations  that  are  coupled  [Eqs. 
(6a)  and  (6b)].  As  a  reference,  we  start  with  the  simulation  of  two  electrons,  initially  at  the 
(0,0)  state  and  the  (0,-1)  state,  without  the  Coulomb  or  exchange  coupling  (Fig.  5).  In  both 
cases,  the  energy  is  evenly  divided  between  kinetic  and  potential.  Figure  6  is  the  same  two 
states,  but  with  only  the  Coulomb  interaction.  (This  corresponds  to  the  case  where  the  particles 
have  opposing  spins.)  The  (0,0)  particle  has  not  changed  substantially  but  the  (0,-1)  has 
effectively  been  pushed  outward,  in  keeping  with  the  Coulomb  repulsion.  This  is  evidenced  by 
the  radius,  which  has  increased,  and  the  fact  that  the  particle  has  substantially  increased  its 
potential  energy  since  it  was  pushed  further  outward  where  the  potential  is  higher.  Note  also 
that,  besides  the  energies  of  the  individual  particles,  there  is  a  Coulomb  energy  due  to  their 
interaction.  Figure  7  shows  the  same  two  particles  with  both  Coulomb  and  exchange 
interactions,  corresponding  to  the  case  where  the  particles  have  the  same  spin.  Notice  that  the 
exchange  interaction  tends  to  pull  the  particles  together  as  evidenced  by  the  increase  in  radius  of 
the  (0,0)  particle  and  decrease  of  the  radius  of  the  (0,-1)  particle  compared  to  Fig.  6.  Note  also 
that  the  Coulomb  energy  has  increased  with  this  increased  closeness,  but  the  exchange  energy  is 
a  negative  quantity  that  reduces  the  total  energy  of  the  system. 

The  recent  development  of  the  single-electron  capacitance  spectroscopy  (SECS) 
technique  allows  experimentalists  to  observe  the  energy  levels  of  individual  electrons  being 
added  to  a  quantum  dot  [3].  This  technique  records  the  energy  needed  to  add  an  electron  under 
the  influence  of  a  magnetic  field,  a  quantity  known  as  the  chemical  potential.  We  will  attempt  to 
use  the  FDTD  method  to  simulate  the  behavior  of  the  first  two  electrons. 

Or  course,  the  first  particle  is  just  the  lowest  level  of  the  Fock-Darwin  level  model  in 
Fig.  4.  It  is  repeated  as  the  dashed  line  labeled  “First  Particle”  in  Fig.  8.  For  the  second 
electron,  we  must  account  for  the  Coulomb  interaction  in  all  cases  and  the  exchange  interaction 
when  the  particles  have  the  same  spins.  For  instance,  we  may  assume  that  the  first  particle  has 
spin  up,  since  this  is  the  lowest  energy  state  when  the  magnetic  field  is  applied.  Therefore,  the 
second  particle  will  start  in  the  (0,0)  state  with  spin  down.  It  will  have  Coulomb  but  no 
exchange  interaction  with  the  first  particle.  The  additional  energy  to  the  system  is  plotted  as  the 
dash-dot  line  labeled  “Singlet”  in  Fig.  8.  Remember,  this  is  the  chemical  potential,  so  it  is  the 
energy  of  the  second  particle  itself,  plus  the  Coulomb  energy.  Suppose  the  second  particle  were 
a  (0,-1)  state  with  spin  up.  It  would  have  a  Coulomb  and  an  exchange  interaction  with  the  first 
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particle.  The  additional  energy  as  a  function  of  the  magnetic  field  is  plotted  as  the  solid  line 
labeled  “Triplet”  in  Fig.  8.  Note  that  these  two  lines  intersect  at  about  LIT,  indicated  by  the 
arrow.  These  two  lines  are  needed  to  simulate  the  second  particle,  because  the  particle  starts  in 
state  (0,0),  spin  down,  and  after  the  magnetic  field  reaches  1.1  T,  it  flips  its  spin  and  goes  to  the 
lower  (0,-1)  state. 

Palecios  et  al.  [5]  carried  out  a  theoretical  study  using  exact  diagonalization  and  the 
unrestricted  Hartree-Fock  approximation  for  up  to  15  electrons.  Their  data,  along  with  that  of 
the  FDTD  simulation,  are  plotted  in  Fig.  9  for  the  first  two  electrons.  (The  two  lines  of  the 
second  particle  from  Fig.  8  have  been  consolidated  to  one,  but  the  arrow  still  marks  the 
crossover  points.)  Clearly,  the  agreement  is  very  good. 

VI.  Conclusion 

We  have  presented  explicit  space-  and  time-domain  Hartree-Fock  simulations  of  two- 
electron  wavepacket  dynamics  in  quantum  dots  based  on  the  FDTD  method.  From  these 
simulations,  we  have  extracted  the  eigenstates  and  eigenenergies.  The  results  are  in  excellent 
agreement  with  the  existing  values  in  the  literature. 

Although  the  calculations  presented  here  do  not  in  and  of  themselves  yield  new  physical 
results,  the  power  of  the  FDTD  technique  must  be  borne  in  mind.  As  an  explicit  space-domain 
technique,  one  can  avoid  difficulties  associated  with  constructing  single-particle  orbitals,  which 
are  used  in  computations  based  on  Slater  determinants  [13].  Although  we  chose  a  parabolic 
confinement  potential,  there  is  no  additional  cost  in  choosing  any  other  potential.  A  Slater- 
determinant-based  calculation  would  require  computing  the  single-particle  orbitals  for  each 
potential  chosen.  Of  course,  the  tradeoff  is  in  the  size  of  the  spatial  mesh  chosen  for  our 
calculations. 

The  calculations  presented  here  are  also  easily  generalized  to  more  particles;  we  estimate 
that  for  the  values  of  the  parameters  used  here,  up  to  four  electrons  should  not  be  too 
computationally  taxing.  In  addition,  simulations  including  spin-orbit  coupling  are  under  way 
and  will  be  the  topic  of  a  future  publication.  Finally,  a  time-domain  approach  allows  one  to 
explore  manifestly  dynamical  properties,  such  as  wavepacket  dynamics  or  the  response  to  time- 
dependent  fields. 
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Figure  Captions 

Figure  1  The  real  part  of  the  wave  functions  of  four  of  the  lowest  eigenstates  of  a  two- 
dimensional  harmonic  oscillator  with  a  ground  state  of  5.4  meV. 

Figure  2  The  time  evolution  of  the  real  part  of  the  wavefunction  of  the  (0,1)  state  through  one 
revival.  Note  that  the  wavefunction  moves  counterclockwise.  The  imaginary  part  of 
this  wavefunction  is  identical,  but  the  lobes  are  90  degrees  out  of  phase  to  the  real 
function  shown  here.  The  wavefunction  of  the  (0,-1)  state  would  move  in  exactly 
the  same  way,  but  in  the  clockwise  directions. 

Figure  3  Contour  plot  of  the  real  part  of  the  waveform  for  the  (0,1)  state  with  an  increasing 
magnetic  field  in  the  perpendicular  direction.  The  increased  B  field  leads  to  a  tighter 
confinement,  as  evidenced  by  the  smaller  radii  and  the  increase  in  energy.  (Energies 
are  in  meV,  radii  are  in  nm.) 

Figure  4  Plot  of  the  Fock-Darwin  levels  for  a  two  dimensional  harmonic  oscillator  with  a 
ground  state  energy  of  5.4  meV. 

Figure  5  The  real  parts  of  the  wavefunctions  of  the  (0,0)  and  (0,-1)  states  of  the  two- 
dimensional  harmonic  oscillator.  On  the  left  side  are  mesh  diagrams  and  on  the  right 
side  are  contour  plots  of  the  same  functions.  These  are  two  independent  particles, 
i.e.,  there  is  neither  a  Coulomb  nor  an  exchange  interactions. 

Figure  6  The  same  two  particles  as  shown  in  Fig.  5,  but  with  the  Coulomb  interactions  added 
(This  corresponds  to  two  particles  with  opposite  spins).  The  Coulomb  repulsion 
has  added  to  the  energy  of  each  particle.  The  radius  of  the  (0,-1)  particle  has 
increased  substantially  because  it  has  been  pushed  outwards.  In  addition  to  the 
kinetic  and  potential  energies  of  each  particle,  there  is  an  energy  of  6.67  meV 
associated  with  the  repulsion  of  the  two  particles  to  each  other. 

Figure?  The  same  two  particles  as  shown  in  Fig.  5,  but  with  both  the  Coulomb  and 
exchange  interactions  added  (The  two  partieles  have  the  same  spins).  The  exchange 
has  decreased  the  energy  of  each  particle.  The  radius  of  the  (0,-1)  particle  has 
moved  a  little  back  towards  the  center.  In  addition  to  the  kinetic  and  potential 
energies  of  each  particle,  there  is  an  energy  of  7.38  meV  associated  with  the 
repulsion  of  the  two  particles  to  each  other  and  one  of  -2.58  meV  associated  with  the 
exchange  force. 
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Figure  8  The  chemical  potential  as  function  of  the  magnetic  field  for  the  first  two  particles  in  a 
quantum  dot.  The  plot  of  the  first  particle  is  simply  the  lowest  level  of  the  Fock- 
Darwin  plot  of  Fig.  4.  The  second  particle,  however,  starts  off  in  the  singlet 
configuration  of  the  (0,0)  state  and  flips  its  spin  to  the  triplet  configuration  to  take 
advantage  of  the  lower  energy  level  of  the  (0,-1)  state  as  the  magnetic  field  increases 
above  LIT. 

Figure  9  Comparison  of  the  simulated  FDTD  data  (as  shown  by  solid  lines)  and  the  calculated 
data  of  Palacious  et  al.  [5]  (as  shown  by  the  asterisks). 
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Figure  2  The  time  evolution  of  the  real  part  of  the  wavefunction  of  the  (0,1)  state 
through  one  revival.  Note  that  the  wavefunction  moves  counterclockwise.  The 
imaginary  part  of  this  wavefunction  is  identical,  but  the  lobes  are  90  degrees  out  of 
phase  to  the  real  function  shown  here.  The  wavefunction  of  the  (0,-1)  state  would 
move  in  exactly  the  same  way,  but  in  the  clockwise  directions. 
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Figure  3.  Contour  plot  of  the  real  part  of  the  waveform  for  the  (0,1)  state  with  an 
increasing  magnetic  field  in  the  perpendicular  direction.  The  increased  B  field  leads 
to  a  tighter  confinement,  as  evidenced  by  the  smaller  radii  and  the  increase  in 
energy.  (Energies  are  in  meV,  radii  are  in  nm.) 


Fock-Darwin  Levels 


Figure  4.  Plot  of  the  Fock-Darwin  levels  for  a  two-dimensional 
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Figure  5.  The  real  parts  of  the  wavefunctions  of  the  (0,0)  and  (0,-1)  states  of  the  two 
dimensional  harmonic  oscillator.  On  the  left  are  the  mesh  diagrams  and  one  the  right 
are  the  contour  plots  of  the  same  functions.  These  are  two  independent  particles,  i.e., 
there  is  neither  a  Coulomb  nor  an  exchange  interactions. 
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Figure  6.  The  same  two  particles  as  shown  in  Fig.  5,  but  with  the  Coulomb  interactions 
added  (This  corresponds  to  two  particles  with  opposite  spins).  The  Coulomb  repulsion 
has  added  to  the  energy  of  each  particle.  The  radius  of  the  (0,-1)  particle  has  increased 
substantially  because  it  has  been  pushed  outwards.  In  addition  to  the  kinetic  and 
potential  energies  of  each  particle,  there  is  an  energy  of  6.67  meV  associated  with  the 
repulsion  of  the  two  particles  to  each  other. 
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Figure  7.  The  same  two  particles  as  shown  in  Fig.  5,  but  with  both  the  Coulomb  and 
exchange  interactions  added  (The  two  particles  have  the  same  spins).  The  exchange 
has  decreased  the  energy  of  each  particle.  The  radius  of  the  (0,-1)  particle  has  moved  a 
little  back  towards  the  center.  In  addition  to  the  kinetic  and  potential  energies  of  each 
particle,  there  is  an  energy  of  7.38  meV  associated  with  the  repulsion  of  the  two 
particles  to  each  other,  and  one  of  -2.58  meV  associated  with  the  exchange  force. 


FDTD  Simulation  of  the  First  Two  Particles 


Figure  8.  The  chemical  potential  as  function  of  the  magnetic  field  for  the  first  two 
particles  in  a  quantum  dot.  The  plot  of  the  first  particle  is  simply  the  lowest  level  of  the 
Fock-Darwin  plot  of  Fig.  4.  The  second  particle,  however,  starts  off  in  the  singlet 
configuration  of  the  (0,0)  state  and  flips  its  spin  to  the  triplet  configuration  to  take 
advantage  of  the  lower  energy  level  of  the  (0,-1)  state  as  the  magnetic  field  increases 
above  1.1  T. 


Chemical  Potentials:  FDTD  vs.  Palacios 


Figure  9.  Comparison  of  the  simulated  FDTD  data  (  as  shown  by  solid  lines)  and 
the  calculated  data  of  Palacious  et  al.  [5]  (as  shown  by  the  asterisks). 
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Abstract 


With  the  present  interest  in  nanostructures,  such  as  quantum  dots,  a  need  exists  to  have 
a  flexible  method  to  be  able  to  determine  eigenvalues  and  eigenstates  for  those  stmctures  that  do 
not  lend  themselves  to  existing  analytical  methods.  In  this  paper  we  present  a  method  that 
accomplishes  this  by  using  a  simulation  of  the  Schroedinger  equation  based  on  the  finite- 
difference  time-domain  (FDTD)  method.  This  method  is  capable  of  simulating  any  structure 
within  the  limits  of  discritization.  By  initializing  a  simulation  with  a  test  function,  the 
eigenfrequencies  are  determined  through  a  Fourier  transform  of  the  resulting  time-domain  data 
collected  at  a  sample  point.  Another  simulation  implements  a  discrete  Fourier  transform  at  the 
eigenfrequencies  at  every  cell  in  the  problem  space,  from  which  the  eigenfunctions  can  be 
constructed. 

Keywords:  Quantum  dots.  Quantum  theory,  FDTD  methods.  Eigenvalues 
and  eigenfunctions.  Semiconductor  device  modeling. 

Classifications:  39-02,  81-02 
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1.  Introduction 


There  is  an  increasing  interest  in  semiconductor  structures  on  the  order  of  hundreds  of 
nanometers  as  potential  computing  devices  and  for  other  applications  [1].  However,  for 
anything  but  the  simplest  canonical  structures,  a  solution  for  the  quantum  mechanical  states  of 
electrons  is  difficult.  For  this  reason,  numerical  rather  than  analytic  approaches  must  be  used  to 
determine  the  eigenenergies  and  eigenfunctions  of  such  structures. 

In  this  paper  we  present  a  formulation  using  the  finite-difference  time-domain  (FDTD) 
method  [2, 3].  This  method  is  widely  used  in  electromagnetic  simulation.  Here  we  present  the 
FDTD  formulation  of  the  Schroedinger  equation;  an  explicit  method  that  begins  with  the  time- 
domain  Schroedinger  equation  and  approximates  the  temporal  and  spatial  derivatives  as 
difference  equations  [3].  In  this  formulation,  the  potential  need  not  be  a  simple  structure  (e.g.  a 
harmonic  oscillator  or  square  well);  any  arbitrary  configuration  can  be  simulated.  Furthermore, 
the  influence  of  a  magnetic  field  is  described.  The  expectation  values  of  observables,  such  as 
kinetic  and  potential  energy,  can  be  easily  calculated. 

The  determination  of  eigenenergies  is  accomplished  by  storing  the  time-domain  data 
from  a  test  point  during  a  simulation  and  determining  the  eigenfrequencies  by  post  processing, 
i.e.,  taking  the  Fourier  transform  of  the  stored  time-domain  data.  The  eigenfunctions  are 
reconstructed  by  rerunning  the  simulation  and  taking  a  discrete  Fourier  transform  at  the 
previously  determined  eigenfrequencies.  This  is  accomplished  by  initializing  the  problem  space 
with  a  test  function  and  then  carrying  out  a  discrete  Fourier  analysis  at  the  eigenfrequencies  at 
every  point  in  the  problem  space  [4,  5].  At  the  end  of  the  simulation,  one  has  the  amplitudes 
and  phases  of  the  eigenfunctions  at  every  point  within  the  problem  space  for  each 
eigenfrequency.  From  this  information,  one  can  determine  the  eigenfunctions. 

We  begin  with  a  description  of  the  FDTD  formulation  of  the  Schroedinger  equation.  As 
an  initial  example,  we  describe  the  simulation  of  a  quantum  dot,  which  is  implemented  by  a 
two-dimensional  harmonic  oscillator  simulation  [6]:  the  third  dimension  can  be  assumed  to  be 
well  confined.  From  the  standpoint  of  the  numerical  technique,  there  is  nothing  special  about  a 
harmonic  oscillator  potential.  We  use  this  to  illustrate  the  method  because  it  has  an  analytic 
solution,  from  which  we  can  verify  the  accuracy  of  the  FDTD  simulation. 

Finally,  we  take  a  “stadium”  potential  to  use  the  full  power  of  the  method  to  determine 
the  low  lying  eigenfunctions  in  a  structure  for  which  there  is  no  analytic  solution. 
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II.  FDTD  Formulation  of  the  Schroedinger  Equation 
Basic  Formulation 

The  time  dependent  Schroedinger  equation  is  given  by  [7] 

VV(^0  +  V(r,t)  •  vr(r,0 

dt  2m 

or 

-jV(r,t)  ■  xif{r,t).  (1) 

dt  2m  n 

To  avoid  using  complex  numbers,  we  will  split  the  variable  into  its  real  and  imaginary 

components: 

V^(r,0  =  WrealM  +  W .  (2) 

Inserting  Eq.  (2)  into  Eq.  (1)  and  separating  the  real  and  imaginary  parts  results  in  the  following 
coupled  set  of  equations: 

^rea^r,t)  ^  (3  a) 

dt  2m  n 

=  (3  b) 

dt  2m  n 

We  begin  by  assuming  a  two-dimensional  space.  We  write  t/r"(i,7)  =  yr(/- Ax,y  •  Ax,n  - At) 
where  i,  j.  and  n  are  indices  and  Ax  and  At  are  the  spatial  and  temporal  steps,  respectively. 
(The  same  spatial  interval  Ax  is  used  in  both  the  x  and  y  directions.)  Starting  with  Eq.  (3a),  the 
finite-difference  approximations  in  space  and  time  result  in 

y/''real(i,j)  -  yr''~\eal{i,j) 

At 

fl  1  y/"  ^'^imag{i  +  l,j)~^y^"  *^^/max(b7)  +  V^"  ^'^imagii  —  I,  j) 

2m  (Axf  _  +  \l/''~^'^imagii,j  +  1)  +  \l/''~^^^imagii  “  1,7  -  1) 

+  —  F(i,  7 )  •  i/r "  inuig  (i,  j) 

n 

from  which  we  get 
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y/  realiijj)  —  realijfj) 

At  h  +  1,7)  “  4y/'”  ^^^imag{i^j)  +  y/  imag(/ ”  Uj) 

Ax^  2m  ^  +  V^”"^'^nu/^(r,7  +  1)  +  ll/''~‘^^^imas(ij  -  1) 

+  — y(/,7)  •  VA” 
n 


The  two  coupled  equations  can  then  be  written 


y/^real(iyj)  —  ^real(i,j) 


At  h  fV^”  ^^^imag 


Ax  2m 


(/  +  1,  j)  —  4y/^  ^^^imag  (i,  7)  +  y/^  imag  {}  1, 7*) 


+  yr'^'^^^^imagiU]  +  1)  +  ^  1) 


n 


j.  _n+l/2 
if/  imag 


0»7)~V^"  imagii,  j) 


At  ti  \^t”real(,i  +  \,j)  —  2y/'^real{i,j)  +  y/^  reali}  l-j) 


Ax^  2m 


+  y/^real  (l,  j  +  l)  +  ^"real  (?,  j  “  1) 


■  ^(*)  j)  ■  V^"re<j/  (j)  7)  • 


Harmonic  Oscillator  Simulation 


One  of  the  well-known  canonical  problems  in  quantum  mechanics  is  the  harmonic 
oscillator.  The  harmonic  oscillator  potential  is  given  by  [7] 

y{x,y)  =  h,\x'^+y^). 

In  the  finite-difference  formulation,  this  is  expressed  as 

v(i,  j)  =  ^ko-[{i-  ic)  ■  Ax +  {j-  jc)  -Axf .  (5) 

It  is  centered  around  (ic,  jc) ,  the  middle  of  the  problem  space.  The  parameter  kg  is 

kQ=m-  (O^ 

and  cOq,  the  frequency  of  oscillation,  is  related  to  the  ground  state  energy  of  the  system  by 

Eo  =  ticoQ.  (6) 
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An  Electron  in  a  Magnetic  Field 


This  section  describes  the  additional  influence  of  a  magnetic  field  oriented  in  the  z 

direction.  We  do  this  by  simulating  the  vector  potential  A  [7] 

A  =  B^i-yi +xj),  (7) 


which  gives 


B  =  VxA  =  1-BqZ 


where  z  is  a  unit  vector  in  the  z  direction. 


The  following  formalism  implements  this  into  the  Schroedinger  equation: 


^  •  a1  \}r{x,y)  +  V{x,y)  •  v^(x,y) 
dt  lm\i  J 

- -qWA  -A-qV  +  q^A^\l/(x,y)  +  V(x,y)  ■  y/ix,y)  • 
.  i  i  J 


(8) 


U 

2m 


The  three  A  terms  used  in  Eq.  (8)  are  determined  by 


A^  =  B,Hy"+x^), 


V7.  (-^ 


Ba{-yi+xj)  =  B„ 


d 


^  3  d 


Putting  these  into  Eq.  (8)  results  in  the  following: 


d\i/ix,y)  _  h 
dt 


2mh 


2  \ 


d_ 

dx^  dx^ 


i\l/(rx,y) 


2mh 

+  V(x,y)-\i/(x,y). 


iy/(x,y)  +  2 


2m 


d  d 

-y - bx — 

dx  dy 


\ir{rx,y) 


(9  a) 
(9  b) 

(9  c) 


Now  we  take  the  finite-difference  approximations  to  the  spatial  and  temporal  derivatives,  along 
with  the  formulation  of  the  harmonic  oscillator  potential  of  Eq.  (5): 
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y/  real{i,j)  —  Ij/  real{i,j) 

h  At  ~4^  imag^i^j) y/  iniag{i yf  imagii~l,j) 

2m  Ax^  _'\‘y/''~^^^imagiij  +  1)  +  y/''~^^^imag(ij  -  1) 

+  ^  AtlAx'^ii  -  icf  +  Ax^U  -  jcf]y/''~''^inmg(i,j) 

2mh  '- 

,  .  Av'‘~^realii,j  +  l)-ll/''''real{iJ-l)) 

Ax{i  -  ic)  ^ - — - 

Br.q  .  2Ax 

H - At 

-  AxO  -  jcy - — - ^ 

I  2Ax  J 

+  *  [(^*  ”  ‘  ^  U  ““  J^)  *  V^”  iniag{i,j)^ 


Of  course,  there  will  be  a  similar  equation  corresponding  to  the  imaginary  part  of  y/ . 


Calculating  the  Observables 

Two  quantities  of  importance  in  quantum  mechanics  are  the  expected  values  of  the 
kinetic  energy  and  the  potential  energy.  They  are  calculated  from  \i/[x,y)  as  follows: 

Kinetic  Energy 

The  expected  value  of  the  kinetic  energy  is  given  by 

The  Laplacian  operator  is  approximated  by 

^  \¥real(i-'^J)-^¥realiiJ)  +  ¥real(^  +  ^J) 

¥reA^,j)  -  _  1)  +  J  +  1) 

The  kinetic  energy  in  the  simulation  is  calculated  by 

{T)  =  '  ’  ¥inu,g(^J)]  '  Vr.a/Cb^)  +  i'^^¥inu.g(U)^  •  (1 1) 

j=l  ,=1 

Potential  Energy 

The  expected  value  of  the  potential  is 

{V)  =  {Y\V\\i/)  =  j_^V(x,  y)\il/{x,  y)f  dr , 

which  is  calculated  by 

(^)  =  E  (12) 

;=1  ,=l 
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III.  Determining  the  Eigenfunctions 

The  first  step  in  determining  the  eigenfunctions  is  to  determine  the  corresponding 
eigenfrequencies.  We  do  this  by  initiating  the  problem  with  a  test  function  at  a  point  in  the 
problem  space  and  storing  the  time-domain  data  at  the  point  of  origin  of  the  test.  The  test 
function  is  chosen  as  a  narrow  pulse  to  insure  that  it  will  be  a  superposition  of  as  many  of  the 
eigenfunctions  as  possible.  After  a  sufficient  number  of  time  steps  (typically  about  30,000  to 
50,000),  the  program  is  halted,  and  a  Fourier  analysis  of  the  time  domain  data  collected  at  the 
source  point  is  performed.  The  peaks  in  the  Fourier  analysis  determine  the  eigenfrequencies. 

Once  we  determine  the  eigenfrequencies,  the  eigenfunctions  can  be  constructed. 
Because  the  Fourier  analysis  resulted  in  a  peak  in  the  Fourier  domain,  we  know  that  the 
eigenfunction  was  present  in  the  original  test  function.  Therefore,  by  performing  another 
Fourier  analysis  at  every  point  in  the  problem  space  at  that  frequency,  the  amplitude  and  phase 
of  that  eigenfunction  can  be  reconstmcted. 

The  Fourier  analysis  to  determine  the  eigenfrequencies  is  accomplished  by  post 
processing,  i.  e.,  an  analysis  of  the  data  takes  place  after  the  simulation  is  complete  at  only  one 
data  point  within  the  problem  space.  To  construct  the  eigenfrequencies,  however,  we  need  to 
perform  the  Fourier  analysis  everywhere  in  the  problem  space,  but  for  a  limited  number  of 
frequencies.  This  is  done  by  a  discrete  Fourier  transform  (DFT)  during  the  simulation, 
sometimes  referred  to  as  a  “running  Fourier  transform”  [4,  5]. 

To  cany  out  the  Fourier  analysis  at  one  frequency,  say  j^,  we  could  compute  the 
following: 


(13) 

0 

(We  take  the  lower  level  of  integration  to  be  zero  because  we  assume  all  functions  are  causal, 
since  we  initialize  them  to  zero  when  the  simulation  begins.)  Taking  Eq.  (13)  into  the  sampled 
time  domain  gives 

N  p  1 

^(n  •  At)  =  Y\Wreai(f^ '  '  [cos(2;5r,n  •  At)  +  ism{27tf^n  ■  Af)] 

n=0 


or 
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N  ^ 

■  At)  =  2,[v're«/("  •  At)co%{2TtfyAt  ■  n)  -  ■  At)  ■  sin(2;»5^Af  ■  n)J  (14  a) 

/i=0 

N 

(N •  At)  =  '^[y/reaiif^ '  At)smi27tf^At •  n)  +  y/i^^{n ■  At)  ■  cos{27if^At ■  n)J.  (14  b) 

Two  equations  like  (14  a)  and  (14  b)  are  needed  for  every  frequency  of  interest.  They 
are  calculated  at  every  cell  in  the  problem  space.  This  is  made  possible  by  noticing  that  the 
summation  in  Eq.  (14a),  for  instance,  can  be  calculated  by  the  following  simple  equation  while 
the  simulation  program  is  running: 

+  (¥„.!("  ■  A0cos(2;i/;A(  ■  n)  -  ■  A()  •  sin(2;5f;  At  ■  n)). 

At  every  point  N,  the  new  value  is  calculated  by  adding  one  more  term  to  the  value  at  N-1. 

When  the  program  has  ran  for  a  sufficient  length  of  time,  the  results  of  Eq.  (14  a)  and 
(14.b)  are  used  to  calculated  the  amplitude  and  phase  at  that  frequency.  These  are  used  to 
reconstruct  the  time-domain  eigenfunctions.  This  process  is  illustrated  in  Fig.  1.  Note  that 
there  is  one  optional  step  indicated  by  the  dashed  line.  Once  the  eigenfunctions  have  been 
determined,  one  of  them  can  be  used  to  initialize  the  FDTD  program.  This  is  useful  to  insure 
that  it  is  indeed  an  eigenfunction  that  remains  stable. 

IV.  Determining  the  Eigenfunctions  of  a  Harmonic  Oscillator  Potential 

In  this  section,  we  will  determine  the  eigenfunctions  and  eigenvalues  of  a  harmonic 
oscillator  potential.  This  is  a  problem  with  an  analytic  solution  which  provides  a  check  of  the 
accuracy  of  the  method.  We  will  use  a  two-dimensional  harmonic  oscillator  potential  with  a 
ground  state  energy  of  2  meV.  The  cells  in  the  simulation  program  are  5  nanometers  (nm)  and 
the  time  steps  are  0.1  femtoseconds  [Ax  and  At,  respectively,  in  Eq.  (10)].  The  problem  space 
is  50x50,  effectively  modeling  an  area  which  is  250x250  nm.  We  begin  by  initializing  the 
problem  with  a  test  function  in  the  middle  of  the  problem  space  (Fig.  2).  This  is  a  Gaussian 
pulse  with  an  envelope  of  15  nm.  Note  the  values  of  kinetic  and  potential  energy,  as  calculated 
by  Eq.  (11)  and  (12).  Since  this  test  particle  is  sitting  on  the  bottom  of  the  potential,  it  initially 
has  very  little  potential  energy.  After  2000  iterations,  or  0.2  picoseconds  (ps),  the  wave 
function  has  begun  to  expand  outward.  It  has  exchanged  some  of  its  kinetic  energy  for  potential 
energy,  but  the  total  energy  remains  the  same.  The  waveform  continues  oscillating  in  and  out. 
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After  50,000  iterations  (5  ps),  the  simulation  is  halted.  Figure  3  shows  the  time  domain  data  at 
the  source  point.  The  Fourier  transform  of  this  data  is  shown  in  Fig.  4a.  The  first  few  peaks 
occur  at  2,  6,  10,  and  14  meV  where  frequencies  have  been  converted  to  energies  by 
multiplication  by  Planck’s  constant. 

We  initiate  the  FDTD  simulation  once  again  with  the  same  test  function.  However,  this 
time  the  program  takes  the  discrete  Fourier  transform  at  the  frequencies  corresponding  to  the 
above-mentioned  energies  (0.48,  1.45,  2.42,  and  3.39  THz).  At  the  end  of  the  simulation,  the 
Fourier  data  is  used  to  reconstruct  the  eigenfunctions  shown  in  Fig.  4  (b). 

The  two-dimensional  harmonic  oscillator  eigenfunctions  can  be  characterized  by  the 
quantum  numbers  n  and  I,  where  n  is  a  positive  integer  corresponding  to  the  number  of  nodes  in 
the  wave  function  as  one  moves  out  radially  from  the  center,  the  principle  quantum  number  [1]. 
The  other  number  /  corresponds  to  angular  momentum,  and  is  clearly  0  for  all  of  these  functions 
since  they  are  all  radially  symmetric.  The  eigenenergies  corresponding  to  the  quantum  numbers 
n,  I  are 

E,,,  =  £„(2/1+|;|+1),  (15) 

where  Eq  is  2  meV  for  the  present  case.  Given  that  I  =  0,  the  four  lowest  states  correspond  to  n 
=  0,  1,  2,  and  3,  which  result  in  energies  of  2,  6,  10,  and  14  meV,  as  we  had  determined  from 
Fig.  4  a. 

That  we  did  not  find  eigenfunctions  for  values  of  /  other  than  zero  should  not  be 
surprising  since  we  started  with  a  symmetric  test  function  in  a  symmetric  potential.  Suppose, 
however,  that  we  begin  with  a  test  function  50  nm  from  the  center  and  initialize  it  with  a  wave 
function  that  has  a  wavelength  of  100  nm  and  a  Gaussian  envelope  of  35  nm  (Fig  5).  Note  that 
this  test  function  starts  with  significant  potential  as  well  as  kinetic  energy.  As  the  simulation 
progresses,  the  function  moves  around  the  potential,  maintaining  constant  values  of  kinetic  and 
potential  energies,  since  it  remains  at  about  a  constant  level  in  the  potential.  After  5  ps,  the 
simulation  is  halted,  and  we  have  the  time-domain  data  at  the  source  point  (Fig.  6).  After  doing 
a  Fourier  transform,  we  have  the  data  of  Fig.  7(a).  The  peak  at  2  meV  will  be  the  same  (n=0, 
1=0)  eigenfunction  of  Fig.  4(b),  so  we  will  run  the  discrete  Fourier  transform  at  4,  6,  8,  and  10 
meV.  The  resulting  eigenfunctions  are  shown  in  Fig.  7(b).  Note  that  these  all  have  n  =  0,  and 
positive  integer  values  of  I  { 2|/|  is  the  number  of  nodes  moving  circumferentially  about  the  dot 
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[1]}.  The  energy  values  correspond  to  those  calculated  by  Eq.  (15)  for  /  =  1,  2,  3,  and  4. 
Obviously,  we  could  run  other  test  functions  that  would  result  in  mixed  values  of  n  and  1. 


Suppose  we  have  the  same  problem  with  a  magnetic  field  of  2  Telsa.  The  magnetic  field 
has  its  greatest  effect  on  states  with  angular  momentum,  so  we  will  rerun  the  previous 
simulation  with  a  particle  of  wavelength  of  100  nm  in  a  Gaussian  envelope  of  35  nm  initialized 
50  nm  from  the  center.  The  resulting  time-domain  function  is  similar  in  form  to  Fig.  6,  but 
results  in  the  Fourier  amplitude  of  Fig.  8(a).  The  pattern  is  the  same,  but  the  peaks  have  been 
shifted.  Specifically,  the  /  =  1  state  has  shifted  from  4  meV  to  6.7  meV,  the  /  =  2  state  from  6  to 
10.9  meV,  etc.  In  fact,  these  are  the  energy  levels  predicted  from  [1]: 


Enj  = 


Ihco,  /  2  -h  /  2f  +  (2n  -h  j/]  -h  l)  , 


(16) 


where  co  ,  is  the  cyclotron  frequency. 


«c=- 


eB 


(17) 


m 


{For  m  in  Eq.  (17),  we  use  the  reduced  mass  of  an  electron  in  GaAs  [1].}  The  resulting 
eigenfunctions  in  Fig.  9(b)  are  similar  to  those  in  Fig.  7(b),  except  they  are  slightly  closer  to  the 
center,  which  results  from  the  additional  confinement  imposed  by  the  B  field. 


V.  Eigenfunctions  of  a  Stadium  Potential 

In  this  section  we  will  use  the  same  method  to  find  eigenfrequencies  and  functions  of  a 
configuration  whose  values  and  functions  are  not  previously  known.  We  will  use  the  “stadium” 
potential  [8]  illustrated  in  Fig.  9.  This  has  a  ground  potential  of  0;  the  stadium  is  shaped  by  a 
potential  barrier  of  1  eV.  As  long  as  the  test  functions  are  on  the  order  of  meV,  we  have  almost 
complete  containment. 

We  begin  by  initializing  a  test  function  at  point  A  and  letting  the  simulation  run  for 
30,000  iterations.  We  then  take  the  Fourier  transform  of  the  time  domain  data  collected  at  the 
source,  which  results  in  the  spectrum  shown  in  Fig.  10a.  From  this,  we  ascertain  that  there  are 
eigenfrequencies  at  0.31,  0.88,  1.79,  and  2.9  meV.  When  we  rerun  the  simulation  computing 
the  discrete  Fourier  transform  at  these  frequencies,  we  obtain  the  functions  of  Fig.  10(b).  Note 
that  we  have  obtained  eigenfunctions  that  are  symmetric  in  both  directions.  This  is  not 
surprising,  since  our  test  function  started  from  the  middle  of  the  stadium. 
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In  order  to  look  for  other  eigenfunctions,  we  repeat  the  simulation  by  initializing  the  test 
function  at  point  B,  which  results  in  the  frequency  response  shown  in  Fig.  11a.  The  lowest 
peak  is  at  0.31  meV;  this  is  obviously  the  ground  state  illustrated  in  Fig.  10,  so  there  is  no  sense 
redoing  this  one.  We  take  the  next  frequencies  which  are  1.05, 1.67,  2.25,  and  2.87  meV.  The 
result  is  the  functions  shown  in  Fig.  1  lb. 


In  order  to  convince  ourselves  that  these  truly  are  eigenfunctions,  we  can  use  them  to 
initialize  the  simulation.  We  will  take  the  eigenfunction  at  2.87  meV.  Using  this,  the  simulation 
was  repeated,  as  illustrated  in  Fig.  12.  As  the  simulation  progresses,  there  is  movement  in  the 
function;  but  at  1.45  picoseconds,  it  has  returned  to  its  original  form.  The  revival  time, 
calculated  by  the  inverse  eigenfrequency 


revival  ^  .00283  CV 


(18) 


predicts  this.  One  can  also  verify  that  is  independent  of  time. 


VI.  Discussion 

We  have  described  and  demonstrated  a  simulation  technique  to  determine  the 
eigenfunctions  for  arbitrary  nanostructures.  The  simulation  is  based  on  the  finite-difference 
time-domain  (FDTD)  method.  The  strength  of  this  method  is  its  flexibility;  there  is  no  limit  to 
the  complexity  of  the  potential  of  the  stmcture  being  simulated,  other  than  what  can  be 
accommodated  by  the  discretization.  This  method  is  also  capable  of  simulating  the  influence  of 
a  magnetic  field.  The  accuracy  of  the  method  was  demonstrated  by  simulating  a  two 
dimensional  harmonic  oscillator  with  and  without  the  magnetic  field.  The  expected  energy 
values  and  eigenfunctions  were  obtained.  The  full  flexibility  of  the  method  was  demonstrated 
by  finding  the  eigenfunctions  of  a  stadium-shaped  potential. 

This  technique  is  powerful  enough  to  study  the  dynamics  in  systems  where  the  time 
dependence  plays  a  central  role,  such  as  for  wavepackets  or  in  the  presence  of  time-dependent 
fields.  In  the  future,  we  will  show  that  the  technique  may  also  be  generalized  in  a 
straightforward  way  to  deal  with  few-electron  wavefunctions  in  quantum  dots  within  the 
Hartree-Fock  approximation. 


12 


Acknowledgements 


This  material  is  based  upon  work  supported  by  the  U.  S.  Army  Research  Office  under  ^ant 
number  DAAH04-96- 1-0406,  and  by  a  grant  for  supercomputer  time  from  the  San  Diego 
Supercomputer  Center.  D.  S.  Citrin  was  supported  by  the  Office  of  Naval  Research  and  by  Ae 
National  Science  Foundation  through  grant  no.  DMR-9705403 


13 


REFERENCES 


1 .  R.  C.  Ashoori,  Electrons  in  artificial  atoms.  Nature,  379,  413  (1996). 

2.  A.  Taflove,  Computational  Electrodynamics— the  Finite-Difference  Time-Domain 
Method.,  (Artech  House,  Boston,  MA,  1995). 

3.  D.  Sullivan,  Electromagnetic  Simulation  Using  the  FDTD  Method.  (IEEE  Press, New 
York  2000). 

4.  C.  M.  Purse,  S.  P.  Mathur,  and  O.  P.  Gandhi,  Improvements  to  the  finite-difference  time- 
domain  method  for  calculating  the  radar  cross  section  of  a  perfectly  conducting  target,  IEEE 

Trans.  Microwave  Theory  and  Tech.,  38,  919  (1990). 

5.  D.  M.  Sullivan,  Mathematical  methods  for  treatment  planning  in  deep  regional 
hyperthermia,  IEEE  Trans.  Microwave  Theory  and  Tech.,  39,  864  (1991). 

6.  J.  Jacak,  P.  Hawrylak,  and  A.  Wojs,  Quantum  Dots,  Springer,  Berlin  1997). 

7.  C.  Cohen-Tannoudji,  B.  Diu,  and  F.  Laloe,  Quantum  Mechanics.  (John  Wiley  and  Sons, 
New  York  (1977). 

8.  S.  Tomsovic  and  E.  J.  Heller,  Long-time  semiclassical  dynamics  of  chaos:  the  stadium 
billiard.  Physics  Review  E,  47,  282  (1993). 


14 


Figure  Captions 

Figure  1  Block  diagram  for  the  procedure  used  to  determine  the  eigenfunctions  of  an  arbitrary 
quantum  nanostructure. 

Figure  2  Time-domain  evolution  of  the  testing  function  initiated  at  the  center  of  the  two- 
dimensional  harmonic  oscillator  potential.  The  original  testing  function  is  a 
Gaussian  envelope  of  3  nm.  Note  that  the  total  energy  of  the  function  remains 
constant,  but  there  is  an  exchange  between  kinetic  and  potential  energy.  (Only  the 
real  part  of  the  wavefunction  is  shown.) 

Figure  3  Time-domain  data  sampled  at  the  source  point  for  the  simulation  illustrated  in  Fig.  2. 

Figured  (a)  Fourier  analysis  of  the  time-domain  data  of  Fig.  3,  (b)  Eigenfunctions 
reconstructed  via  the  DFT  at  the  eigenfrequencies  found  in  (a). 

Figure  5  Time-domain  evolution  of  the  testing  function  initiated  50  nm  from  the  center  of  the 
harmonic  oscillator  potential.  This  function  was  initiated  as  a  waveform  with  a 
wavelength  of  100  nm  inside  a  Gaussian  envelope  of  35  nm,  giving  it  substantial 
kinetic  energy.  (Only  the  real  part  of  the  wavefunction  is  shown.) 

Figure  6  Time-domain  data  sampled  at  the  source  point  for  the  simulation  illustrated  in  Fig.  5. 

Figure?  (a)  Fourier  analysis  of  the  time-domain  data  of  Fig.  6,  (b)  Eigenfunctions 

reconstmcted  via  the  DFT  at  the  eigenfrequencies  found  in  (a). 

Figure  8  (a)  Fourier  analysis  of  the  time-domain  data  collected  for  a  test  function  moving  in  a 

magnetic  field  of  2  Telsa.  Note  that  the  form  of  the  Fourier  analysis  is  very  similar 
to  Fig.  7a,  but  the  peaks  are  shifted,  (b)  Eigenfunctions  reconstructed  via  the  DFT 
at  the  eigenfrequencies  found  in  (a). 

Figure  9  Diagram  of  a  “stadium”  potential.  The  inside  of  the  stadium  is  at  zero  potential, 
while  the  outside  is  1  eV. 

Figure  10  (a)  Fourier  analysis  of  the  time-domain  data  at  source  point  A  for  the  potential 

illustrated  in  Fig.  9,  (b)  Eigenfunctions  reconstructed  via  the  DFT  at  the 

eigenfrequencies  found  in  (a). 

Figure  1 1  (a)  Fourier  analysis  of  the  time-domain  data  at  source  point  B  of  Fig.  9  for  a  testing 
function  initiated  at  source  point  B.  (b)  Eigenfunctions  reconstructed  via  the  DFT  at 
the  eigenfrequencies  found  in  (a). 

Figure  12  Time  evolution  of  a  function  in  the  stadium  potential  that  was  initiated  from  the  2.87 
meV  eigenfunction  of  Fig.  11b. 


Figure  1.  Block  diagram  of  the  procedure  used  to  determine 
the  eigenfunctions  of  an  arbitrary  quantum  nanostructure. 
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Figure  2.  Time-domain  evolution  of  the  testing  function  initiated  at  the  center  of  the  two- 
dimensional  harmonic  oscillator  potential.  The  original  testing  function  is  a  Gaussian  pulse  of 
15  nm.  Note  that  the  total  energy  of  the  function  remains  constant,  but  there  is  an  exchange 
between  kinetic  and  potential  energy.  (Only  the  real  part  of  the  wavefunction  is  shown.) 
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Figure  3.  Time-domain  data  sampled  at  the  source  point  for  the  simulation  illustrated 
in  Fig.  2. 
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Figure  5.  Time-domain  evolution  of  the  testing  function  initiated  50  nm  from  the 
center  of  the  harmonic  oscillator  potential.  This  function  was  initiated  as  a  waveform 
with  a  wavelength  of  100  nm  inside  a  Gaussian  envelope  of  35  nm,  giving  it 
substantial  kinetic  energy.  (Only  the  real  part  of  the  wavefunction  is  shown.) 


Figure  6.  Time-domain  data  sampled  at  the  source  point  for  the  simulation  illustrated 
in  Fig.  5. 
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Figure  8.  (a)  Fourier  analysis  of  the  time  domain  collected  for  a  test  function  moving  in 
a  harmonic  oscillator  with  a  magnetic  field  of  2  T.  Note  that  the  form  of  the  Fourier 
analysis  is  very  similar  to  Fig.  7a,  but  the  peaks  are  shifted,  (b)  Eigenfunction 
reconstmcted  via  the  DFT  at  the  eigenfrequencies  found  in  (a). 
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Figure  9.  Diagram  of  a  “stadium”  potential.  The  inside  of  the  stadium  is  at  zero 
potential,  while  the  outside  is  at  1  eV. 
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Figure  10.  (a)  Fourier  analysis  of  the  time-domain  data  at  source  point  A  for  the  potential 
illustrated  in  Fig.  9;  (b)  Eigenfunctions  reconstmcted  via  the  DFT  at  the  eigenfrequencies  found 
in  (a). 
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Figure  12.  Time  evolution  of  a  function  in  the  stadium  potential  that  was  initiated  from  the  2.87 
meV  eigenfunction  of  Fig.  lib. 
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Abstract 

The  need  for  short,  high-frequency  pulses  for  spectoscopy  and  imaging  has  motivated 
the  search  for  methods  to  shape  these  pulses.  Conductive  apertures  are  often  used  for 
spatiotemporal  shaping  of  THz  pulses.  Computer  simulation  can  facilitate  the  design  of  such 
apertures  by  allowing  one  to  evaluate  different  configurations  and  their  related  parameters. 
Three-dimensional  simulations  using  the  FDTD  method  with  techniques  described  here  make 
this  possible  without  the  need  for  extraordinary  computer  resources.  The  accuracy  of  this 
method  is  confirmed  by  comparison  with  previously  published  experimental  data.  The 
versatility  of  the  method  is  demonstrated  by  the  simulation  of  apertures  of  various  shapes  and 
sizes. 
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Supercomputer  Center.  D.  S.  Citrin  was  supported  by  the  Office  of  Naval  Research  and  by  the 
National  Science  Foundation  through  grant  no.  DMR-9705403. 


1.  Introduction 


Over  the  past  few  years,  terahertz  (10”'^  5"')  technology  has  been  experiencing 
explosive  growth.  The  growth  in  this  area  has  been  fueled  largely  by  the  need  for  faster  signal 
processing  and  communications,  for  high-resolution  spectroscopy,  atmospheric  and 
astrophysical  remote  sensing,  and  for  imaging.  A  number  of  methods  for  the  generation  of  THz 
radiation  have  been  developed.  Many  of  these  methods  use  a  photoconductive  antenna  in  which 
an  ultrashort  optical  pulse  is  incident  upon  a  strongly  biased  semiconducting  material  [1].  Other 
methods  use  heterostructures,  such  as  p-i-n  junctions  [2]  and  strained  superlattices  [3].  Some 
of  the  applications  using  these  radiated  THz  pulses  require  a  pulse  with  a  specific  shape  to 
improve  the  prospects  for  optical  control  over  material  behavior  and  structure.  A  method  for 
spatiotemporal  shaping  of  the  THz  pulses  by  a  diffraction  through  conductive  apertures  in  a 
two-dimensional  space  has  recently  been  reported  [4].  Simulation  of  the  THz  pulse  shaping  is 
desirable  because  it  gives  an  opportunity  to  find  the  kind  of  pulse  shape  that  can  be  obtained 
using  different  apertures.  Simulation  also  allows  the  investigation  of  the  influence  of  variations 
in  the  parameters  much  faster  and  more  efficiently  than  experimental  methods. 

This  paper  investigates  the  effects  of  the  apertures  in  a  three-dimensional  problem  space. 
Three-dimensional  simulation  is  more  desirable  than  two-dimensional  simulation  for  three 
primary  reasons.  First,  the  shape  in  the  third  dimension  can  be  varied.  Second,  it  is  desirable 
to  observe  the  simulated  pulse  anywhere  in  the  three-dimensional  space  and  not  just  be  limited 
to  the  single  horizontal  plane  provided  by  a  two-dimensional  simulation.  Third,  the  wave 
physics  of  a  two-dimensional  space  is  different  from  that  of  one-  and  three-  dimensions  [5].  In 
one-  and  three-dimensional  waves,  the  shape  of  the  wave  is  preserved  as  it  travels  from  its 
initial  position,  whereas  in  a  two-dimensional  wave  there  exists  a  “wake”  that  trails  off  after  the 
pulse  as  it  travels  forward  in  time. 

To  simulate  the  far  field  of  the  aperture  in  a  three-dimensional  problem  space, 
prohibitively  large  amounts  of  core  computer  memory  would  normally  be  required.  To 
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overcome  this,  a  transformation  is  used  to  determine  the  far  field  from  the  near  field  in  the 
aperture.  Furthermore,  the  problem  space  can  be  reduced  by  exploiting  symmetry. 

2.  Simulation  Methods 


2.1  The  Finite-Difference  Time-Domain  Method 

The  finite-difference  time-domain  (FDTD)  method  [6]-[8]  is  one  of  the  most  popular 
numerical  methods  for  solving  problems  in  electromagnetics.  The  time-domain  Maxwell’s 
equations  in  free  space  are  given  by 


dE 

dt 


=  -^Vx£r 

Co 


(la) 


dH 

dt 


--VxE. 


The  FDTD  formulation  implements  these  as  difference  equations,  two  of  which  are 


(lb) 


/1+1/2 


+  - 


At-Cn 


Ax 


H/ii  + 1  /  2,7, A:  + 1  /  2)  -  H/(i  -1/ 2,j,k  +  l  /  2) 

+ 1  /  2,fc  + 1  /  2)  +  - 1  /  2,A:  + 1  /  2)J 


(2  a) 


+ 1  /  2,7  + 1  /  2,k)  =  H;(i  +  l/2,j  +  l/2,k) 


At -Cn 


Ax 


'E;ii +ij +1/ 2,  k)-  E;(i,j +\  /  2,  k) 
-E;  {i  +  \l2,j  +  \,k)  +  E;{i  +  \l2, 7,  k) 


(2  b) 


Similar  equations  follow  for  E^,  E^,  H^,md  Hy  fields. 

Figure  1(a)  is  a  diagram  of  a  typical  experimental  setup  for  spatiotemporal  filtering  of 
THz  pulses,  and  Fig  1(b)  shows  the  region  to  be  simulated.  The  problem  domain  must  be 
broken  into  cells  for  simulation  in  the  FDTD  method.  (We  use  cubes  with  the  same  dimension 
in  each  direction.)  Some  thought  must  be  given  to  the  size  of  the  cells.  We  will  be  dealing  with 
pulses  that  have  frequency  components  as  high  as  1  THz.  This  corresponds  to  a  wavelength  of 
300  |J,m.  To  maintain  accuracy  under  the  simulation  conditions  that  we  are  using,  it  is  desirable 
to  have  a  sampling  rate  of  twenty  points  per  wavelength,  i.  e.,  a  cell  size  of  15  jim.  The  high 
sampling  rate  reduces  the  propagation  error  [6]  and  also  reduces  error  due  to  the  “staircasing” 
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effect  that  occurs  when  modeling  metal  objects  with  edges  that  do  not  lie  in  one  of  the 
rectangular  directions.  This  means  that  to  simulate  a  very  small  area  of  4  mm  by  7  mm  by  4 
mm,  a  problem  space  of  200x500x200  cells  is  required.  Implementing  an  FDTD  program  of 
this  size  would  require  400  megawords  of  core  memory  in  the  computer. 

To  overcome  this  problem,  two  techniques  can  be  used:  The  symmetry  of  the  problem 
can  be  utilized  to  reduce  the  simulated  space  to  one  fourth,  and  a  transformation  can  be  used 
whereby  the  far  field  can  be  determined  by  the  fields  in  the  aperture,  eliminating  the  need  to 
simulate  the  far  field. 

2.2  Symmetry 

We  use  the  fourfold  symmetry  of  the  aperture  to  reduce  the  computational  problem 
space.  The  four  quadrants  are  symmetric  about  the  x  =  0  and  y  =  0  planes;  therefore,  only  one 
of  the  four  quadrants  is  simulated.  This  technique  has  been  used  before  in  FDTD  simulation 

[9]. 

2.3  The  Near-Field  to  Far-Field  Transformation 

The  diffracted  pulse  from  the  aperture  is  observed  in  the  radiating  far  field.  To  include 
the  far  field  in  our  computation,  a  large  amount  of  computer  memory  is  required;  however,  by 
implementing  the  near-field  to  far-field  transformation  described  in  a  previous  paper  [10],  the 
required  computational  space  is  dramatically  reduced.  In  this  transformation,  the  FDTD- 
computed  time  domain  data  at  the  end  of  the  aperture  is  used  to  calculated  the  far-field  values 
(Fig.  2).  A  brief  description  follows. 

We  start  with  the  vector  potential 


(3) 


where  Jf  =  |r  -  r| .  The  H  field  due  to  the  electric  current  density  J  is  obtained  by 
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H  =  VxA 


1 

=  — V  X  f7(r’)^ — dV 
4;r  ^  i 


where  r  =  r^x  +  r^y  +  =  j^i  1^1  =  -  Jc'  )^  +  (y  -  y'f  +  z^f^-  Using  the  duality  principle, 

the  E  field  due  to  the  magnetic  current  density  M  is  given  by 


£  =  — jJ(M(r')xr)  ;^  + 


RJ  R 


For  an  x-polarized  incident  plane  wave,  the  direction  normal  to  the  aperture  is  the  z  direction,  so 


M  =  2E.xn 


■  2E^xxz=--2Ej, 


where  E^  is  the  field  in  the  aperture.  Hence, 


ys  /s  /  /V  y^\ 

yxr  =  [r^z-r^x). 


Therefore,  the  x-polarized  E  field  is  given  by 


* = ^  iJ  V^' . 


Taking  this  into  the  sampled  time  domain,  we  get  the  following  equation: 


KTl 


E:{i,j\n-2-RJ 


where  EJ'it  ,f  ,n-2-Rj^)  is  the  time-retarded  field  at  the  end  of  the  aperture,  and 
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del_E:{,i  ,f  ,n-2RJ  =  E;{i 


(11) 


More  details  are  given  in  [10]. 

Since  symmetry  is  being  used,  for  every  point  in  the  plane  at  the  end  of  the  aperture, 
there  exist  three  other  points  in  the  other  three  quadrants  that  must  be  considered  to  obtain  the 
true  transformation.  As  can  be  seen  from  Fig.  3,  the  point  P,,  which  is  in  the  simulated 
quadrant,  has  symmetry  points  Pj,  P3,  and  P4.  The  E  field  values  in  the  aperture  at  the  four 
points  are  the  same,  but  the  distances  to  the  far-field  observation  point  to  each  of  these  four 
symmetric  points  is  not  the  same.  Hence,  when  implementing  the  transformation  along  with  the 
symmetry,  this  has  to  be  taken  into  account. 

Using  the  two  techniques  described  above,  we  reduce  the  sample  problem  that  required 
a  domain  of  200x400x200  cells  to  a  domain  of  100x200x100  cells  that  requires  only  26 
megawords  of  core  computer  memory.  The  simulated  problem  domain  is  shown  in  Fig.  4.  The 
near-field  to  far-field  transformation  reduces  the  problem  space  to  a  small  area  around  the 
aperture,  and  the  use  of  symmetry  further  reduces  the  problem  space  by  seventy-five  percent. 
Note  that  every  side  of  the  problem  space,  except  the  symmetry  planes,  has  a  twelve  cell 
perfectly  matched  layer  (PML)  to  absorb  outgoing  scattered  waves  [11]. 

3.  Verification 

To  verify  the  accuracy  of  the  results  obtained  from  the  three-dimensional  simulation 
using  the  symmetry  and  the  transformation  techniques,  we  compare  the  resulting  pulses  from 
the  full  FDTD  simulation  with  those  obtained  using  the  transformation  techniques.  A  Gaussian 
pulse  with  a  spread  of  0.375  picoseconds  (or  15  time  steps)  is  diffracted  through  a  rectangular 
aperture  of  thickness  1.00  mm  and  width  of  0.5  mm.  Comparisons  are  made  at  selected  points 
[Fig.  1(b)].  Figures  5(a)  and  5(b)  show  the  comparison  on-axis  at  distances  of  0.24  mm  and 
1 .5  mm  respectively  from  the  aperture.  Figures  6(a)  and  6(b)  show  the  comparison  at  off-axis 
distances  of  0.45  mm  in  the  y-direction  and  x-direction,  respectively. 
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To  verify  the  method  against  experimental  data,  we  simulate  the  setup  described  by 
Bromage,  et  al  [4],  who  generated  an  incident  pulse  by  the  illumination  of  a  GaAs  wafer  with  a 
Ti:sapphire  laser  operating  at  810  nm.  The  resulting  pulse  is  shown  in  Fig.  7(a).  This  was 
filtered  by  a  rectangular  aperture  of  thickness  1.7  mm  and  a  width  of  0.5  mm.  For  comparison, 
the  resulting  pulse  and  the  filtering  aperture  were  simulated.  The  results,  shown  in  Fig.  7(b), 
are  in  good  agreement.  The  reconstructed  incident  pulse  in  Fig.  7(a)  was  used  for  all 
subsequent  simulations  in  this  paper. 

4.  Aperture  Simulation 

In  this  section,  we  use  the  three-dimensional  simulation  program  to  evaluate  the  pulse 
shaping  characteristic  from  different  apertures.  We  begin  with  a  rectangular  aperture  and 
evaluate  the  effects  of  the  various  parameters.  We  then  look  at  other  tapered  and  circular 
apertures  for  comparison. 

The  specific  parameters  of  the  rectangular  aperture  are  the  thickness  (1)  in  the  z  direction, 
the  width  (d)  in  the  y  direction,  and  the  height  (h)  in  the  x  direction  (Fig.  4).  All  the  observation 
points  in  this  section  are  at  a  distance  of  7  mm  from  the  aperture.  The  results  for  the  pulses 
obtained  by  varying  the  thickness  and  the  width  of  the  aperture  were  found  to  be  in  accordance 
with  the  results  shown  in  the  paper  by  Bromage,  et  al.  [4].  When  d  =  0.5  mm  and  h  =  0.15 
mm,  we  observe  that  as  the  thickness  increased,  considerable  pulse  shaping  occurs  (Fig.  8). 
Furthermore,  the  ringing  increases  due  to  the  fact  that  the  aperture  acts  as  a  resonating  cavity 
with  increasing  thickness.  Also  for  I  of  0.23  mm  and  /i  of  0.15  mm,  as  the  width  increases, 
narrower  pulses  are  obtained  as  shown  in  Fig.  9.  This  is  can  be  explained  by  the  following: 
for  a  rectangular  waveguide  the  width  d,  the  cut-off  frequency  in  free  space  is  given  by  [12]: 


As  the  height  of  the  aperture  is  increased,  we  found  that  only  the  amplitude  of  the  resulting 
pulse  increased. 


7 


The  three-dimensional  capability  was  particularly  useful  in  simulating  shapes  that  could 
not  adequately  be  represented  in  two-dimensions.  Now  a  tapered  aperture  can  be  simulated,  as 
shown  in  Fig.  10.  The  resulting  pulse  after  passing  through  a  tapered  aperture  like  Fig.  10  is 
shown  in  Fig.  11(a).  Another  simulation  used  an  aperture  tapered  on  both  edges;  the  resulting 
pulse  is  shown  in  Fig.  11(b).  The  results  of  the  tapered  apertures  show  pulses  that  die  out 
smoothly  without  excessive  ringing,  as  opposed  to  the  untapered  rectangular  apertures. 

The  use  of  a  circular  aperture  was  also  investigated.  After  a  process  of  varying 
parameters,  similar  to  the  procedure  using  the  rectangular  aperture,  we  found  that  a  conical 
shaped  aperture  (Fig.  12)  gave  the  pulse  shapes  of  Fig.  13.  These  shapes  are  very  appealing, 
particularly  Fig.  13(a),  because  they  present  a  couple  of  very  smooth  cycles  without  ringing. 
However,  they  do  nothing  to  eliminate  the  secondary  pulse. 

Which  of  the  above  pulses  is  best  depends  upon  the  application.  It  was  our  intention  to 
demonstrate  the  versatility  of  the  method  and  the  wide  variety  of  pulses  that  can  be  obtained 
with  various  apertures. 

5.  Conclusion 

In  this  paper,  we  carried  out  three-dimensional  simulations  of  THz  spatiotemporal  pulse 
shaping.  The  computations  exploit  the  symmetry  of  the  problem  and  a  transformation  that  lets 
one  determine  the  far  field  from  the  fields  in  the  aperture.  We  also  demonstrated  that  three- 
dimensional  simulation  provides  a  more  complete  description  than  two-dimensional  simulation 
because  three-dimensional  simulation  allows  the  simulation  of  all  parameters  in  a  wider  variety 
of  structures.  Such  simulations  facilitate  the  exploration  of  pulse  shapes  achievable  in  complex, 
though  technologically  feasible,  geometries.  In  particular,  the  full  vector  nature  of  the 
simulations  allows  for  the  accurate  presentation  of  pulse  shapes  far  from  the  propagation  axis 
where  scalar  techniques  are  known  to  fail. 
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Captions 

Filtering  of  THz  pulses:  (a)  Experimental  setup,  (b)  Problem  domain  to  be 
simulated.  The  parameters  I  and  d  are  the  length  and  thickness  of  the  aperture, 
respectively. 

The  electric  field  in  the  aperture  E{x' ,y')  can  be  used  to  calculate  the  far  field 
Ey{x,y).  The  aperture  A  is  an  opening  in  a  metal  plate  that  lies  in  the  XY  plane. 

The  distances  from  the  four  symmetric  points  in  the  plane  containing  the  aperture  to 
the  far-field  observation  point  are  not  the  same,  even  though  the  E  field  values  at  the 
four  points  are  equal. 

Diagram  of  the  reduced  FDTD  problem  space.  By  using  the  near-  to  far-field 
transformation,  any  point  in  the  far  field  can  be  calculated  by  the  values  in  the 
aperture.  Furthermore,  by  using  symmetry,  only  one-fourth  of  the  problem  is 
simulated. 

Comparison  of  the  waveforms  determined  by  the  FDTD  simulation  directly  (solid 
line)  versus  those  determined  by  the  transformation  (dashed  line).  A  plane  wave 
Gaussian  pulse  of  0.375  picoseconds  was  passed  through  a  rectangular  aperture 
with  a  thickness  of  1  mm  and  a  width  of  0.5  mm.  The  waveforms  are  at  observation 
points  perpendicular  to  the  center  axis  (z  direction)  at  distances  of  0.24  mm  (a)  and 
1.5  mm  (b). 

Comparison  of  the  waveforms  determined  by  the  FDTD  simulation  directly  (solid 
line)  versus  those  determined  by  the  transformation  (dashed  line).  A  plane  wave 
Gaussian  pulse  of  0.375  picoseconds  was  passed  through  a  rectangular  aperture 
with  a  thickness  of  1  mm  and  a  width  of  0.5  mm.  The  waveforms  are  at  observation 
points  1.5  mm  away  in  the  z  direction,  but  offset  4.5  mm  in  the  y-direction  (a)  and 
x-direction  (b). 

Comparison  of  simulated  versus  the  experimental  data  of  Bromage,  et  al.  [4].  (a)  is 
the  incident  pulse  and  (b)  is  the  resulting  waveform  after  a  rectangular  aperture  of 
1.7  mm  thickness  and  0.5  mm  width. 

Resulting  waveforms  obtained  by  varying  the  thickness  of  a  rectangular  aperture 
with  a  width  of  0.5  mm  and  a  height  of  0.15  mm. 

Resulting  waveforms  obtained  by  varying  the  gap  width  of  a  rectangular  aperture 
with  a  thickness  of  0.225  mm  and  a  height  of  0. 1 5  mm. 
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Figure  10.  A  tapered  aperture. 

Figure  11.  Waveforms  resulting  from  apertures  tapered  on  one  side  (a)  similar  to  Fig.  10,  or 
tapered  on  both  sides  (b). 

Figure  12.  A  circular  aperture. 

Figure  13  Waveforms  resulting  from  circular  apertures:  (a)  major  radius  0.225  mm  and  a 
thickness  of  0.15  mm,  (b)  major  radius  0 .525  mm  and  thickness  of  0.3  mm. 
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Figure  1.  Filtering  of  THz  pulses:  (a)  Experimental  set  up,  (b)  Problem  domain  to  be  simulated. 
The  paramters  I  and  d  are  the  length  and  thickness  of  the  aperture,  respectively. 


Figure  2.  The  electric  field  in  the  aperture  E(x’,y')  can  be  used  to 
calculate  the  far  field  Ef(X,y,z).  The  aperture  A  is  an  opening  in 

a  metal  plate  that  lies  in  the  XY  plane. 


Figure  3.  The  distances  from  the  four  symmetric  points  in  the  plane 
containing  the  aperture  to  the  far-field  observation  point  are  not  the 
same,  evan  though  the  E  field  values  at  the  four  points  are  equal. 
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FDTD  program  in  the  aperture 
are  used  to  determine  the  far 
field  values. 

Figure  4.  Diagram  of  the  reduced  FDTD  problem  space.  By  using 
the  near-  to  far-field  transformation,  any  point  in  the  far  field  can  be 
calculated  by  the  values  in  the  aperture.  Furthermore,  by  using 
symmetry,  only  one-forth  of  the  problem  is  simulated. 
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Figure  5.  Comparison  of  the  waveforms  determined  by  the  FDTD  simulation  directly 
(solid  lines)  versus  those  determined  by  the  transformation  (dashed  line).  A  plane  wave 
Gaussian  pulse  of  0.375  picoseconds  was  passed  through  a  rectangular  aperture  with  a 
thickness  of  1  mm  and  a  width  of  0.5  mm.  The  waveforms  are  at  observation  points 
perpendicular  to  the  center  axis  (z  direction)  at  distances  of  0.24  mm  (a)  and  1.5  mm  (b). 
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Figure  6.  Comparison  of  the  waveforms  determined  by  the  FDTD  simulation  directly 
(solid  lines)  versus  those  determined  by  the  transformation  (dashed  line).  A  plane  wave 
Gaussian  pulse  of  0.375  picoseconds  was  passed  through  a  rectangular  aperture  with  a 
thiclcness  of  1  mm  and  a  width  of  0.5  mm.  The  waveforms  are  at  observation  points  1.5 
mm  away  in  the  z  direction,  but  offset  .45  mm  in  the  y  direction  (a)  and  x-direction  (b). 
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Figure  7.  Comparison  of  simulated  versus  experimental  data  of  Bromage,  et  al.  [4].  (a)  is  the 
incident  pulse  and  (b)  is  the  resulting  wvaeform  after  a  rectangular  aperture  of  1.7  nom  thickness 
and  0.5  mm  width. 
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Figure  8.  Resulting  waveforms  obtained  by  varying  the  thickness  of  a  rectangular  aperture 
witha  width  of  0.5  mm  and  a  heightof  0.15  mm. 
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Figure  9.  Resultingwaveforms  obtained  by  varying  the  gap  width  of  a  rectangular  aperture  with 
a  thickness  of  0.225  mm  and  a  height  of  0. 15  mm. 
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Figure  1 1 .  Waveforms  resulting  from  aperuters  tapered  on  one  side  (5 
tapered  on  both  sides  (b). 
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Figure  12.  A  circular  aperture. 


V/m  V/m 


#  .  ► 


Conical  Aperture 


(a) 


Figure  13.  Waveforms  resulting  from  circular  apertuers:  (a)  major  radius  2.25  mm  and 
a  thickness  of  0.15  mm,  (b)  major  radius  5.25  mm  and  thickness  of  0.3  mm. 


